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1. Introduction 

Simulating strongly interacting fermions continues to be a major challenge in computational 
physics. The standard procedure to deal with fermionic degrees of freedom is to integrate out the 
fermionic fields in order to obtain the fermion determinant det D, where D denotes the Dirac op- 
erator. However, this procedure is not unproblematic. Consider for example a fermion interacting 
with a bosonic field U. After integrating out the fermion fields one obtains det D(U) which yields 
an effective action non-local in the bosonic field. The standard method is now to re-express the de- 
terminant using bosonic 'pseudo-fermions' and use the Hybrid Monte Carlo algorithm [jl|] which in 
essence encodes the non-locality of the fermion determinant in the inverse D(U)~ l . Another prob- 
lem is that the standard approach suffers from critical slowing down (CSD) towards the chiral limit. 
In that limit the correlation length of the fermionic two-point function diverges. As a consequence 
the Dirac operator D(U) develops very small modes and eventually the inverse D{U)~ l becomes 
ill-conditioned. Yet another problem concerns the phase of detD which for Wilson fermions is in 
general non-zero. Hence a probabilistic interpretation of the integration measure, necessary for any 
Monte Carlo simulation, is not possible and leads to a sign problem when an odd number of Wilson 
fermion flavours is simulated. 

Here we propose a novel approach [Q] circumventing the above mentioned problems. It is 
based on the exact hopping expansion of the fermion action, i.e. a reformulation of the fermion 
system as a statistical closed loop model. We develop a simulation algorithm which samples di- 
rectly the fermionic two-point function and in this way eliminates CSD. Moreover, it allows to 
specify the fermionic boundary conditions a posteriori, i.e. after the simulation, and allows simula- 
tions directly in the massless limit. The approach is applicable to the Gross-Neveu (GN) model in 
D = 2 dimensions, to the Schwinger model in the strong coupling limit in D = 2 and D = 3 dimen- 
sions, to supersymmetric quantum mechanics and the N = I and 2 supersymmetric Wess-Zumino 
model in D = 2 dimensions. In the present proceedings we concentrate on the application to the 
GN and the Schwinger model. 

Finally, we would like to emphasise that the reformulation based on the hopping expansion is 
not new [|l [|, ||, ^]. Mostly, however, it has been applied to staggered fermions in the strong cou- 
pling limit where a reformulation in terms of monomers and dimers ^ allows efficient algorithms 



[g, |9p that were subsequently applied to many interesting systems [ |10[ , [TT| , 12], see also the recent 
review by Chandrasekharan [|T|]. For Wilson fermions on the other hand the loop formulation 
has been developed for the Schwinger model in the strong coupling limit and the GN model 



, 15] and what we propose in [Q] is just a very efficient algorithm for these loop formulations. 



2. Loop formulation of Wilson fermions 

We start with the reformulation of D = 2 fermionic systems involving Wilson fermions in terms 
of a statistical loop gas model. We use the GN model, a prototype for strongly interacting fermions, 
as an illustrative example. The model is most naturally formulated in terms of Majorana fermions. 
Employing the Wilson lattice discretisation for a Majorana fermion the Euclidean Lagrange density 
reads 

& = \?ny& - \d*d+*n)Z ~ J ($ T n) 2 (2-1) 
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where £ is a real 2-component Grassmann field, = — ^ r is the charge conjugation matrix and 
d,d* and d are the forward, backward and the symmetric lattice derivative, respectively. In the 
continuum, the massless model enjoys a discrete chiral symmetry £ — > 75^ which on the lattice 
is broken explicitly by the Wilson term ^d*d. The symmetry can be restored in the continuum 
by fine tuning m — » ra c . Further we note that a pair of Majorana fermions may be considered 
as one Dirac fermion, i.e. y = 1/a/2(^i + it^), W = 1 /V2(^[ — i^J)^, exposing the 0(2N) 
flavour symmetry explicitly. Since integrating out Majorana fermions yields the Pfaffian of the 
antisymmetric Dirac operator, the model with 2N Majorana fermions is equivalent to N Dirac 
fermions through the identity (PfD) 2N = (detD) N . 

At non- vanishing coupling g ^ one usually employs a Hubbard-Stratonovich transformation 
and introduces the scalar field a oc t, Tc ^t,. WithM(x) = 2 + m + a(x) andP(±jU) = j(lTY(i) the 
action then becomes the sum of monomer and hopping terms 

Scx = ~£e(x)VM(x)S(x)-£e(x)VP(LiMx + fi). (2.2) 
Using the nil-potency of Grassmann elements one can now expand the Boltzmann factor and per- 



form an exact hopping expansion for the Majorana Wilson fermions [|15[]. We emphasise that this 
can be done for any fermionic theory (bilinear in the fermionic fields). At each site, the fields Z, Tc € 
and t, must be exactly paired in order to give a non- vanishing contribution to the path integral, 

/ 9% n( M w/ 2 ^w^w) wW n(^ r w^(M)^(^+A)) i ' ,i{j:) (2.3) 

J x x,p. 

where the occupation numbers m(x) = 0, 1 for monomers and bn[x) = 0, 1 for bonds (or dimers) 
satisfy the constraint 

mto + iEM*) = 1 - ( 2 - 4 ) 
1 n 

This constraint encodes that only closed, non-intersecting paths survive the integration and we end 
up with a closed loop representation of the partition function in terms of monomers and dimers, 
i.e. Z = Y,e C0(£). The weight CO of each loop I can be calculated analytically [Q, ^, 15, |16|] yield- 



ing \oo{i)\ = 2~ c l 2 where c is the number of corners in the loop, while the phase of oo(£) de- 
pends on the geometrical shape of £ In D = 2 dimensions and for a torus geometry of the lattice, 
sign[ft)(^)] G { — 1, 1} depends on the boundary conditions (BC) G {0, 1} and on the number 
of loop windings in direction pL, 

sign[co(£)] = (-1)%( £ m+%) . (2.5) 

As a consequence the overall sign of a given configuration depends only on the fermionic BC and 
the total winding number / = {/^} (modulo 2). 

If we separate all configurations into the equivalence classes J^ ; - where the subscripts i,j 
specify the total winding numbers (modulo 2) in the two directions, then the partition function 
summing over all non-oriented, self-avoiding loops with positive weight, 

Z= £ l«Min M M> ^€JSfooUJSfioUJS%iUJSfn, (2.6) 
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Figure 1: N — 1 Majorana GN model on a 128 2 lattice. Left: Comparison of simulation results (symbols) 
and analytic calculations (dashed lines) for the partition function ratios /Z. The inset shows the repro- 
duction of the zero mode of at m c = 0. Right: Integrated autocorrelation time of the condensate at the 
critical point m c = fitted by %a ~ L z with z = 0.31 (4). The inset shows a fit to a logarithmic dependence 
on L. 



represents a system with unspecified fermionic BC while systems with specific fermionic BC can 
be constructed a posteriori by taking the signs of each class according to 

Z\ = 1Z^- £(-l)^"Z^. (2.7) 

Finally we note that if one considers N > I Majorana flavours the occupation numbers m,^ are 
decorated by the flavour index a and one considers N different loop flavours. The monomer weight 
M(x) depends on the local fermion density £ a m a (x) only and one ends up with a model of locally 
coupled loops. 

In the Schwinger model the hopping term contains a U(l) phase coming from the gauge field 
0u(x), and the non-oriented (Majorana) bonds carry an additional factor °< cosh(^(jc)). More- 
over the gauge field introduces an interaction between the two Majorana flavours proportional to 
±sinh(0„ (x)), These additional factors introduce a sign problem since each loop can now have an 
arbitrary sign. However, in the strong coupling limit, the two flavours are bound together. In the 
present formulation it means that two different Majorana loops lay on top of each other and the 
resulting double loop describes the world line of the bosonic bound state. It also turns out that 
all the signs cancel in a non-trivial way and so the bosonisation is realised explicitly. Eventually 



we end up with a model of non-oriented loops [ 14 ] in which all the loop and monomer weights 



are squared compared to the GN model. Note further that eq.(2.7) no longer applies because the 
fermionic BC have no impact on the BC of the corresponding bosonic bound state - instead the 
relevant partition function is the one where all topological classes contribute positively, i.e. Z. 

3. Simulation algorithm for loops and strings 

A standard procedure to simulate loop gas models as the one described above is to perform 



local loop updates involving plaquette moves only [17, 18]. One problem with such an algorithm 
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is that it can not change between the topological classes J5Joo,-£fio,-£oi,JS?n. Moreover, if the 
correlation length of the system grows large these algorithms become highly inefficient and suffer 



from CSD. Our proposal [g] (subsequently worked out in []19[]) follows the one of Prokof 'ev and 
Svistunov [|(]] and enlarges the configuration space by open fermionic strings. In the GN model 
an open string corresponds to the insertion of a Majorana fermion pair {£, (x) , ^ T (y)^} at position 
x and y into the path integral, and the open string samples directly the correlation function 

G(x,y) = J ne- SG ^(x)^(y) T ^. (3.1) 

This is the reason why CSD is eliminated: configurations are updated on all length scales up to 
0(Q where £ is the correlation length corresponding to the fermionic two point function. As 
a consequence the update remains efficient even at a critical point where the correlation length 
diverges. Contact with the partition functions Z^,. is made each time the open string closes and this 
provides the proper normalisation for the expectation value of the 2-pt. function, (x)% {y) Tc &)z = 
G(x,y)/Z, or any other observables. In practice, the ends of the open string are updated with a 
standard local Metropolis or heat bath procedure [0]. Similar ideas have been around for a long 



time in various other contexts [EQ, 21 , 22] - what is new here is the practical application to Wilson 



fermions and the demonstration that CSD is essentially eliminated. 



4. Absence of critical slowing down 

Before investigating the efficiency of the algorithm, we demonstrate its correctness by compar- 
ing simulation results with analytically know expressions. For this purpose we use the N =1 Majo- 
rana GN model. This model is essentially a free fermion model and can be solved exactly by calcu- 
lating Pfaffians in momentum space. In the left plot of Figure [j] we show the results for the partition 
function ratios Z_j*../Z on a 128 2 lattice from 2M closed path configurations (symbols) compared to 
the exact results (dashed lines). The inset shows the combination = Z^ m — Z_jf 10 — Zg> 01 —Z& u 
which has a zero mode at the critical point m c = 0. The algorithm is indeed able to reproduce the 
zero mode without problems. In order to investigate the efficiency of the algorithm at the criti- 
cal point we measure the condensate (£ r< «f£)zg. The right plot of Figure [I] shows the integrated 
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(m - m } L 



Figure 3: The Schwinger model in the strong coupling limit. Left: Finite size scaling of Z^°/Z for a second 
order phase transition in the universality class of the Ising model. Right: Integrated autocorrelation time 
of the condensate at the critical point m c fitted by Ta ~ L z with z = 0.25(2). The inset shows a fit to a 
logarithmic dependence on L. 

autocorrelation time %a of the condensate as a function of the linear system size L. The dynamic ex- 
ponent z relevant for CSD, i.e. x A ~ L z , turns out to be z — 0.31(4). A dependence logarithmically 
on L can also be fitted to L > 32 yielding -14.2(2.5) +7.1(6) ln(L) with x 2 /dof = 0.18. 

Next we consider the Schwinger model in the strong coupling limit g — > oo as a non-trivial 
example for strongly interacting fermions. In the left plot of Figure |2] we show the partition function 
ratio Zf/Z on various lattices up to L = 512. As in the Majorana GN model we find a zero 
of the partition function which depends only very little on the extent of the lattice. We can use 
Z| (m c ) = as a definition for the critical point m c . It can be determined by a linear fit and 
we obtain m c = —0.686506(27) (cf. right plot in Figure |J) from our simulations on the largest 
lattice with L = 512. Further improvement could be achieved by employing standard reweighting 
techniques as done in [ |23| ] where they obtained m c = —0.6859(4). These calculations indicated a 
second order phase transition in the universality class of the Ising model (with critical exponent 
V ~ 1). Our results in the left plot of Figure [3] now confirm this by demonstrating that the partition 
function ratios Z^°/Z as a function of the rescaled mass (m — m c )L v with v = 1 beautifully collapse 
onto a universal scaling curve. The efficiency of the algorithm and the fact that CSD is essentially 
absent is demonstrated in the right plot of Fig. || where we show the integrated autocorrelation time 
%a of the energy as a function of the linear system size L at the critical point m = m c . The functional 
dependence on L can be well fitted {% 2 /&of = 1.28) by x A ~ L z all the way down to our smallest 
system size L = 8. We obtain z = 0.25(2) which is consistent with just using the largest two system 
sizes. The autocorrelation time may also depend logarithmically on L and a fit to L > 32 yields 
— 13.8(1.9) + 6.6(4) ln(L) with # 2 /dof = 1.00. In any case it is an amazing result that our local 
Metropolis-type update appears to have a dynamical critical exponent close to zero. 

5. Conclusions 

In conclusion, we have presented a new type of algorithm for Wilson fermions in two di- 
mensions. It relies on sampling directly 2-point correlation functions and essentially eliminates 
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critical slowing down. We have successfully tested our algorithm on the Majorana GN model and 
on the Schwinger model in the strong coupling limit and found remarkably small dynamical criti- 
cal exponents. The algorithm definitely opens the way to simulate efficiently generic loop models 
(with positive weights) in arbitrary dimensions, in particular the GN model with any number of 
flavours, the Thirring model, the Schwinger model and QED3 in the strong coupling limit, as well 
as fermionic models with Yukawa-type scalar interactions like the N = 1 and 2 Wess-Zumino mod- 
els, all with Wilson fermions. 
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